home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Programming / Source / Lyapunov / lyap.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-07-31  |  9.5 KB  |  443 lines
  1. /*
  2.  * by Garrett Wollman based on a function posted to Usenet by Ed Kubaitis
  3.  */
  4.  
  5. /* This program is Copyright (C) 1991, Garrett A. Wollman.  This
  6.  * program may be modified and distributed for any purpose and without
  7.  * fee, provided that this notice remains on all such copies
  8.  * unaltered.  Binary distributions not including this source file are
  9.  * prohibited.  Modified versions must be marked with the name of the
  10.  * modifier and the date of modification.
  11.  *
  12.  * Under no circumstances shall Garrett A. Wollman or the University
  13.  * of Vermont and State Agricultural College be held liable for any
  14.  * damages resulting from the use or misuse of this program, whether
  15.  * the author is aware of such possibility or not.  This program is
  16.  * warranted solely to occupy disk space.
  17.  *
  18.  * I'm sorry for all this legal garbage, but it is sadly necessary
  19.  * these days.
  20.  */
  21.  
  22. #include <math.h>
  23. #include <time.h>
  24. #include <stdlib.h>
  25. #include <stdio.h>
  26.  
  27. /* other multiprocessing hosts please add your info here */
  28.  
  29.  
  30. /*
  31.  * Global variables
  32.  */
  33. #ifdef MULTIPROC
  34. int nprocs = 1;
  35. #endif
  36.  
  37. char *whoami;
  38.  
  39. extern int getopt(int,char **,const char *);
  40. extern char *optarg;
  41. extern int optind;
  42. extern int opterr;
  43. extern void perror(const char *);
  44.  
  45. #ifdef SUN_BROKEN_STDLIB
  46. extern volatile void exit(int);
  47. extern void *malloc(int);
  48. extern int atoi(const char *);
  49. #endif
  50.  
  51. extern int strlen(const char *);
  52.  
  53. #ifdef NO_STDIO_PROTOS
  54. extern int fwrite(const char *,int,int,FILE *);
  55. extern int fputc(int,FILE *);
  56. extern int fputs(const char *,FILE *);
  57. extern int pclose(FILE *);
  58. extern void fclose(FILE *);
  59. extern int fprintf(FILE *,const char *,...);
  60. extern int sscanf(char *,const char *,...);
  61. #endif
  62.  
  63. #ifdef MULTIPROC
  64. #define OPTIONS "r:c:v:m:d:s:g:5pt#:o:l:"
  65. #else
  66. #define OPTIONS "r:c:v:m:d:s:g:5pto:l:"
  67. #endif
  68.  
  69. int rows = 512,cols = 512;
  70. double amin,amax;
  71. double bmin,bmax;
  72. double aincr,curA;
  73. double bincr;
  74.  
  75. #define nColors 256        /* don't even try to change this */
  76. struct {
  77.   unsigned char red;
  78.   unsigned char green;
  79.   unsigned char blue;
  80. } colormap[nColors];
  81.  
  82. int Settle = 600;
  83. int Dwell = 2000;
  84. int *Rvec;
  85. double minLya = -5;
  86. int showpos = 0;
  87. double initGuess = 0.5;
  88.  
  89. /*
  90.  * Colormap functions
  91.  */
  92.  
  93. /*
  94.  * set up our initial shades of grey
  95.  */
  96. void init_colormap(void) {
  97.   int i;
  98.   for(i=0; i < nColors; i++)
  99.     colormap[i].red = colormap[i].green = colormap[i].blue = i;
  100. }
  101.  
  102. void read_colormap(const char *name) {
  103.   FILE *cmap;
  104.   char buf[255];        /* magic number */
  105.   int a,b,c,i;
  106.   
  107.   cmap = fopen(name,"r");
  108.   if(!cmap) {
  109.     perror(name);
  110.     exit(-1);
  111.   }
  112.  
  113.   i = 0;
  114.   while(!feof(cmap) && i < nColors) {
  115.     fgets(buf,sizeof buf,cmap);
  116.     /*
  117.      * unparseable lines are comments
  118.      * as is anything after the <r g b> value on a single line
  119.      *
  120.      * per Fractint 15.1 standard
  121.      */
  122.     if(sscanf(buf,"%d %d %d",&a,&b,&c) != 3)
  123.       continue;
  124.  
  125.     colormap[i].red = a;
  126.     colormap[i].green = b;
  127.     colormap[i].blue = c;
  128.     i++;
  129.   }
  130.   fclose(cmap);
  131. }
  132.  
  133.  
  134. /*
  135.  * Lyapunov exponent and coloring calculation
  136.  *
  137.  * original function by Ed Kubaitis, used by permission
  138.  */
  139.  
  140. /*
  141.  * Speedup...
  142.  *
  143.  * Erick B. Wong noticed that $\log_2 d(x_1)+\dots +\log_2 d(x_n)$
  144.  * is, by a principle from high-school algebra, the same as
  145.  * $\log_2 (d(x_1) + \dots + d(x_n))$.  We take advantage of this
  146.  * by unrolling the dwell loop by four (hence the rounding in main())
  147.  * and accumulating by multiplication before taking the log.  Since
  148.  * log() is extremely expensive, this should lead to a considerable
  149.  * speedup, while still allowing complete flexibility in selecting
  150.  * dwell values.
  151.  *
  152.  * Presumably, one might unroll this loop even further, with a 
  153.  * concomitant increase in textual complexity.
  154.  */
  155.  
  156. int lya(double a,double b) {
  157.   double t = 0;
  158.   double r, lya;
  159.   int d;
  160.   double acc;
  161.   double x;
  162.  
  163.   x = initGuess;        /* initialize x */
  164.  
  165.   for(d=0; d < Settle; d++) {
  166.     r = (Rvec[d]) ? b : a;
  167.     x = r*x*(1-x);
  168.   }
  169.  
  170.   for(d=0; d <= Dwell; d+= 4) {
  171.     r = (Rvec[d]) ? b : a;
  172.     x = r * x * (1-x);
  173.     acc = fabs(r - 2*r*x);
  174.  
  175.     r = (Rvec[d+1]) ? b : a;
  176.     x = r * x * (1-x);
  177.     acc *= fabs(r - 2*r*x);
  178.  
  179.     r = (Rvec[d+2]) ? b : a;
  180.     x = r * x * (1-x);
  181.     acc *= fabs(r - 2*r*x);
  182.  
  183.     r = (Rvec[d+3]) ? b : a;
  184.     x = r * x * (1-x);
  185.     t += log(acc * fabs(r - 2*r*x));
  186.  
  187.     if(fabs(x-0.5) < 1e-20) {
  188.       d += 3;
  189.       break;
  190.     }
  191.   }
  192.   lya = (t * 1.442695041)/d;
  193.  
  194.   if(showpos)            
  195.     lya *= -1;
  196.   if(lya < minLya)        /* sanity check! */
  197.     lya = minLya;
  198.  
  199.   if(lya < 0) {
  200.     return (int)(lya / minLya * (double)nColors);
  201.   } else {
  202.     return 0;
  203.   }
  204. }
  205.  
  206.  
  207. /*
  208.  * Usage message
  209.  */
  210. volatile
  211. void usage(void) {
  212.   fprintf(stderr,
  213.       "%s: usage:\n\n"
  214.       "%s [-r rows] [-c cols] [-v program] [-p] [-t]\n\t"
  215.       "[-o file] [-l colormap]\n\t"
  216.       "[-m minLya] [-d Dwell] [-s Settle] [-g initGuess]\n\t"
  217.       "[-5] amin amax bmin bmax bits\n",whoami,whoami);
  218.   exit(-1);
  219. }
  220.  
  221. int main(int argc,char **argv) {
  222.   int c,i,j;
  223.   char *viewprogram = NULL;
  224.   char *outname = NULL;
  225.   FILE *outfile;
  226.  
  227.   char *bits;
  228.  
  229.   double curB;
  230.   int onRow,onCol;
  231.   int *oneRow;
  232.  
  233.   int usepgm = 0;
  234.   int outtext = 0;
  235.  
  236.   init_colormap();
  237.   opterr = 0;
  238.  
  239.   whoami = *argv;
  240.   while((c = getopt(argc,argv,OPTIONS)) != EOF) {
  241.     switch(c) {
  242.     case 'r':
  243.       rows = atoi(optarg);
  244.       if(!rows)
  245.     rows = 512;
  246. /*
  247.     usage();
  248. */
  249.       break;
  250.  
  251.     case 'c':
  252.       cols = atoi(optarg);
  253.       if(!cols)
  254.     cols = 512;
  255. /*
  256.     usage();
  257. */
  258.       break;
  259.  
  260.  
  261.     case 'o':
  262.       outname = optarg;
  263.       viewprogram = NULL;
  264.       break;
  265.  
  266.     case 'v':
  267.       viewprogram = optarg;
  268.       outname = NULL;
  269.       break;
  270.  
  271.     case 'l':
  272.       read_colormap(optarg);
  273.       break;
  274.  
  275.     case 'm':
  276.       if(!sscanf(optarg,"%lf",&minLya))
  277.     usage();
  278.       break;
  279.  
  280.     case 'd':
  281.       Dwell = atoi(optarg);
  282.       if(!Dwell)
  283.     usage();
  284.       /*
  285.        * By forcing the Dwell to be a multiple of four, we can use the
  286.        * speedup noted in lya().  We *could* just do this silently,
  287.        * by if it might make a difference to someone...
  288.        */
  289.       if(Dwell % 4) {
  290.     fprintf(stderr,"Dwell value of %d being rounded off to %d.\n",
  291.         Dwell,
  292.         Dwell += 4 - (Dwell % 4));
  293.       }
  294.       break;
  295.  
  296.     case 's':
  297.       Settle = atoi(optarg);
  298.       if(!Settle)        /* for settle close to 0, use 1 */
  299.     usage();
  300.       break;
  301.  
  302.     case 'g':
  303.       if(!sscanf(optarg,"%lf",&initGuess))
  304.     usage();
  305.       break;
  306.  
  307.     case '5':
  308.       usepgm = 1;
  309.       break;
  310.  
  311.     case 'p':
  312.       showpos = 1;
  313.       break;
  314.  
  315.     case 't':
  316.       outtext = 1;
  317.       break;
  318.  
  319.     default:            /* includes '?' (and '#' when !MULTIPROC) */
  320.       usage();
  321.     }
  322.   }
  323.  
  324.   if(!argv[optind])
  325.     usage();
  326.  
  327.   /*
  328.    * Note to language nit-pickers... The code below is legal because
  329.    * the || operator is a sequence point.
  330.    */
  331.   if(!sscanf(argv[optind++],"%lf",&amin) || !argv[optind] ||
  332.      !sscanf(argv[optind++],"%lf",&amax) || !argv[optind] ||
  333.      !sscanf(argv[optind++],"%lf",&bmin) || !argv[optind] ||
  334.      !sscanf(argv[optind++],"%lf",&bmax) || !argv[optind])
  335.     usage();
  336.  
  337.   bits = argv[optind];
  338.  
  339.   if(Dwell < Settle)
  340.     Dwell = Settle;        /* a bit of sanity */
  341.  
  342.   j = strlen(bits);
  343.  
  344.   Rvec = (int *)malloc((Dwell+1) * sizeof *Rvec);
  345.   if(!Rvec) {
  346.     fputs("Cannot allocate space for Markus vector!\n",stderr);
  347.     exit(-1);
  348.   }
  349.  
  350.   for(i=0; i <= Dwell; i++) {
  351.     if(bits[i % j] == 'a')
  352.       Rvec[i] = 0;
  353.     else if(bits[i % j] == 'b')
  354.       Rvec[i] = 1;
  355.     else {
  356.       fprintf(stderr,"Invalid character \\%o ('%c') in bit string\n",
  357.           (int)(unsigned char)bits[i % j], bits[i % j]);
  358.       exit(-1);
  359.     }
  360.   }
  361.  
  362.   oneRow = (int *)malloc(cols * sizeof *oneRow);
  363.   if(!oneRow) {
  364.     fputs("Cannot allocate space for single-row buffer!\n",stderr);
  365.     exit(-1);
  366.   }
  367.  
  368.   aincr = (amax - amin) / rows;
  369.   bincr = (bmax - bmin) / cols;
  370.  
  371.   if(viewprogram) {
  372.     outfile = popen(viewprogram,"w");
  373.     if(!outfile) {
  374.       fprintf(stderr,"Couldn't popen(\"%s\",\"w\")!\n",viewprogram);
  375.       exit(-1);
  376.     }
  377.   } else if(outname) {
  378.     outfile = fopen(outname,"w");
  379.     if(!outfile) {
  380.       perror(outname);
  381.       exit(-1);
  382.     }
  383.   } else {
  384.     outfile = stdout;
  385.   }
  386.  
  387.   /*
  388.    * Print P[PG]M header...
  389.    */
  390.  /* fprintf(outfile,"P%d %d %d\n",(usepgm?2:3) + (outtext?0:3),cols,rows);
  391.   fprintf(outfile,"%d\n",nColors-1);*/
  392.   putc('\0',outfile);  putc('\0',outfile);
  393.   putc('\0',outfile);  putc('\0',outfile);
  394.   putc('\0',outfile);  putc('\0',outfile);
  395.   putc('\0',outfile);  putc('\0',outfile);
  396.   putc('\0',outfile);  putc('\0',outfile);
  397.   putc('\0',outfile);  putc('\0',outfile);
  398.   fprintf(stderr, "Sizes: %d %d %d %d\n",cols%256,cols/256,rows%256,rows/256);
  399.   putc((unsigned char)(cols%256),outfile);
  400.   putc((unsigned char)(cols/256),outfile);
  401.   putc((unsigned char)(rows%256),outfile);
  402.   putc((unsigned char)(rows/256),outfile);
  403.   putc('\0',outfile);  putc('\0',outfile);
  404.  
  405.  
  406.   for(onRow = 0, curA = amax; onRow < rows; curA -= aincr, onRow++) {
  407.     for(onCol = cols - 1, curB = bmax; onCol; curB -= bincr, onCol--) {
  408.       oneRow[onCol] = lya(curA,curB);
  409.     }
  410.     
  411.     for(i=0; i < cols; i++) {
  412.       unsigned char c;
  413.       c = (unsigned char)oneRow[i];
  414.       if(outtext) {
  415.     if(usepgm)
  416.       fprintf(outfile,"%d ",c);
  417.     else
  418.       fprintf(outfile,"%d %d %d ",colormap[c].red,colormap[c].green,
  419.           colormap[c].blue);
  420.       } else if(usepgm) {
  421.     fwrite((char *)&c,sizeof c,1,outfile);
  422.       } else {
  423.       putc(colormap[c].red,outfile);
  424.       putc(colormap[c].green,outfile);
  425.       putc(colormap[c].blue,outfile);
  426.     /*fwrite((char *)&colormap[c],sizeof colormap[0],1,outfile);*/
  427.       }
  428.     }
  429.     if(!(onRow % 16)) {
  430.       if(outtext)
  431.     fputc('\n',outfile);
  432.       fprintf(stderr,"Done row %d/%d \n",onRow,rows);
  433.     }
  434.   }
  435.  
  436.   if(viewprogram)
  437.     pclose(outfile);
  438.   else if(outname)
  439.     fclose(outfile);
  440.  
  441.   return 0;
  442. }
  443.